###########################################################################################
###### creating shape file versions of outputs ############################################
###########################################################################################

#load packages:
library(R.utils)
library(raster)
library(maptools)
library(rgdal)
library(sp)
library(sf)

#Main working dir:
MyDir<-paste("/home/jc217070/Maxent",sep="")  #define working directory on HPC
OutDirFin<-file.path(MyDir,"Outputs3_Final")
OccDir<-file.path(MyDir,"Maxent_SWDs/spp_occs")
IndSppOutDir<-paste(OutDirFin, "/Final_Individual_Spp", sep="")

#set project name:
projectName<-"WHOSnakes"

#define selection criteria of interest for further analyses:
crit<-c("permutation.importance")   #c("permutation.importance","contribution","Test.gain.without","Test.gain.with.only","AUC.without","AUC.with.only")

#read species reference list:
DatSum<-read.table(file=file.path(paste(MyDir, "/Maxent_LUTables/",projectName,".csv",sep="")), header=TRUE, sep=",") #load data summary for all species
spp<-DatSum$gen_sp_subtax[DatSum$Include=="YES"]
spp<-c("Echis_carinatus","Daboia_russelii","Bungarus_caeruleus","Naja_naja")

#outcats<-c("Current","Median","Q90","Q10")
outcats<-c("Current")

#################################################################################
#make temporary temp dir for rasters in this job and set from default to this dir
dir.create(paste(MyDir,"/RastoShape",sep=""),showWarnings=FALSE)
rasterOptions(tmpdir=paste(MyDir,"/RastoShape",sep=""))
#################################################################################

for (cat in outcats){
  
  ###################################
  # general mapspp info #############
  ###################################
  
  for (sp in spp){
    
    n<-which(DatSum$gen_sp_subtax==sp) #gets row number for current species
    m<-which(spp==sp) #gets row number for current species
    mapsp<-DatSum$ModUnit[n]
      
      
      
      if(file.exists(paste(OccDir,"/",sp,"_allRecs.csv",sep=""))){
      sprecso<-read.table(paste(OccDir,"/",sp,"_allRecs.csv",sep=""), sep=",", header=TRUE)   #read in species recs
      sprecs<-sprecso[sprecso$uncertainty<= (100000),]
      if(length(sprecs[,1])==0){
        sprecs<-sprecso
      }
      
      if(length(sprecs[,1])>0){     #only continue if there's any points
        coords<-data.frame(long=sprecs$long,lat=sprecs$lat)  # get coordinates of occurrence records
        spcoords <- SpatialPoints(coords[,1:2])              # makes spatial version of occurrence coordinates
        spcoords@proj4string<-CRS("+init=epsg:4326")         # define projection to match my .asc files
        
        Sys.sleep(0.1)
        print(paste(Sys.time(), "creating shapefiles for species", m, "out of", length(spp),":", sp, sep=" "))
        
        Sys.sleep(0.1)
        print(paste(Sys.time(), "1. reading original rasters", sep=" "))
        
        cutSDM<-raster(paste(IndSppOutDir,"/",sp,"_",cat,"_CropThreshCDCut.asc",sep=""), crs="+init=epsg:4326")
        cutSDMBin<-raster(paste(IndSppOutDir,"/",sp,"_",cat,"_CropThreshCDCutBin.asc",sep=""), crs="+init=epsg:4326")
        lin_dist.ras<-raster(paste(IndSppOutDir,"/",sp,"_",cat,"_RelCostDistance.asc",sep=""), crs="+init=epsg:4326")
        snapraster<-cutSDM
        snapraster[!is.na(snapraster)] <- NA

        #####################################################################################
        #Known areas of occurrence (buffers around occurrences cut to suitable habitat)
        
        Sys.sleep(0.1)
        print(paste(Sys.time(), "2a. creating 'known to occur' areas", sep=" "))
        
        Buf<-buffer(spcoords, width=25000, dissolve=TRUE) #make buffer around current habitat
        BufKTO<-rasterize(x=Buf,y=snapraster,as.numeric(1))
        KTO<-BufKTO * cutSDM
        values(KTO)<-ifelse(values(KTO)<0.0000001, NA, values(KTO))
        if(length(values(KTO)[!is.na(values(KTO))])==0){
          KTO<-BufKTO
          values(KTO)<-ifelse(values(KTO)<0.0000001, NA, values(KTO))
        }
        KTOBin<-KTO
        values(KTOBin)<-ifelse(values(KTOBin)<0.0000001, NA, 5)
        
        Sys.sleep(0.1)
        print(paste(Sys.time(), "2b. coverting to shapefile", sep=" "))
        
        KTOshp<-rasterToPolygons(KTOBin, dissolve=TRUE)
        KTOshp<-st_as_sf(KTOshp)
        names(KTOshp)[1]<-"VALUE"
        KTOshp$TYPE<-rep("Known Areas of Occurrence",nrow(KTOshp))
        KTOshp$SPECIES<-rep(gsub("_"," ",sp),nrow(KTOshp))
        KTOshp$MODUNIT<-rep(gsub("_"," ",mapsp),nrow(KTOshp))
        KTOshp$CODE_SSN<-rep(DatSum$SSN[n],nrow(KTOshp))
        KTOshp$DATETIME<-rep(date(),nrow(KTOshp))
        KTOshp$CAT<-rep(cat,nrow(KTOshp))
        st_write(KTOshp, paste(IndSppOutDir,"/",sp,"_",cat,"_KTO.shp",sep=""),delete_layer=TRUE)
        
        #########################################
        #CD quantiles:

        Sys.sleep(0.1)
        print(paste(Sys.time(), "3a. creating Cost Distance quantiles", sep=" "))

        lin_dist.ras<-lin_dist.ras*cutSDMBin
        myvals<-values(lin_dist.ras)
        myvals<-myvals[myvals>0]
        CDQ75<-round(quantile(myvals,0.75, na.rm=TRUE),2)
        CDQ75<-ifelse(CDQ75==1,0.99,CDQ75)
        CDQ75<-ifelse(CDQ75< 0.03,0.03,CDQ75)
        print(paste("CDQ75 =",CDQ75))
        CDQ50<-round(quantile(myvals,0.5, na.rm=TRUE),2)
        CDQ50<-ifelse(CDQ50==1,CDQ50-0.02,CDQ50)
        CDQ50<-ifelse(CDQ50==CDQ75,CDQ50-0.01,CDQ50)
        CDQ50<-ifelse(CDQ50< 0.02,0.02,CDQ50)
        print(paste("CDQ50 =",CDQ50))
        CDQ25<-round(quantile(myvals,0.25, na.rm=TRUE),2)
        CDQ25<-ifelse(CDQ25==1,CDQ25-0.03,CDQ25)
        CDQ25<-ifelse(CDQ25==CDQ75,CDQ25-0.02,CDQ25)
        CDQ25<-ifelse(CDQ25==CDQ50,CDQ25-0.01,CDQ25)
        CDQ25<-ifelse(CDQ25<0.0000001,0.0000001,CDQ25)
        print(paste("CDQ25 =",CDQ25))
        CDqs<-lin_dist.ras

        if(is.na(CDQ75)){
          CDQ75<-0.75
          CDQ50<-0.5
          CDQ25<-0.25
          CDqs<-KTO
          print(paste("default CDQ breaks = 0, 0.25, 0.5, 0.75, 1",CDQ75))
        }

        values(CDqs)<-ifelse(values(CDqs)>Q75, 4, values(CDqs))
        values(CDqs)<-ifelse(values(CDqs)>Q50 & values(CDqs) <= Q75, 3, values(CDqs))
        values(CDqs)<-ifelse(values(CDqs)>Q25 & values(CDqs) <= Q50, 2, values(CDqs))
        values(CDqs)<-ifelse(values(CDqs)>0.0000001 & values(CDqs) <= Q25, 1, values(CDqs))
        values(CDqs)<-ifelse(values(CDqs)<=0.0000001, NA, values(CDqs))

        Sys.sleep(0.1)
        print(paste(Sys.time(), "3b. coverting to shapefile", sep=" "))

        CDqsshp<-rasterToPolygons(CDqs, dissolve=TRUE)
        CDqsshp<-st_as_sf(CDqsshp)
        names(CDqsshp)[1]<-"VALUE"
        CDqsshp$TYPE<-rep("Cost Distance Quantile Categories",nrow(CDqsshp))
        CDqsshp$SPECIES<-rep(gsub("_"," ",sp),nrow(CDqsshp))
        CDqsshp$MODUNIT<-rep(gsub("_"," ",mapsp),nrow(CDqsshp))
        CDqsshp$CODE_SSN<-rep(DatSum$SSN[n],nrow(CDqsshp))
        CDqsshp$DATETIME<-rep(date(),nrow(CDqsshp))
        CDqsshp$CAT<-rep(cat,nrow(CDqsshp))
        st_write(CDqsshp, paste(IndSppOutDir,"/",sp,"_",cat,"_CDQs.shp",sep=""),delete_layer=TRUE)

        #########################################
        #suitability quantiles:

        Sys.sleep(0.1)
        print(paste(Sys.time(), "4a. creating SDM quantiles", sep=" "))

        myvals<-values(cutSDM)
        myvals<-myvals[myvals>0]
        SDMQ75<-round(quantile(myvals,0.75, na.rm=TRUE),2)
        SDMQ75<-ifelse(SDMQ75==1,0.99,SDMQ75)
        SDMQ75<-ifelse(SDMQ75< 0.03,0.03,SDMQ75)
        print(paste("SDMQ75 =",SDMQ75))
        SDMQ50<-round(quantile(myvals,0.5, na.rm=TRUE),2)
        SDMQ50<-ifelse(SDMQ50==1,SDMQ50-0.02,SDMQ50)
        SDMQ50<-ifelse(SDMQ50==SDMQ75,SDMQ50-0.01,SDMQ50)
        SDMQ50<-ifelse(SDMQ50< 0.02,0.02,SDMQ50)
        print(paste("SDMQ50 =",SDMQ50))
        SDMQ25<-round(quantile(myvals,0.25, na.rm=TRUE),2)
        SDMQ25<-ifelse(SDMQ25==1,SDMQ25-0.03,SDMQ25)
        SDMQ25<-ifelse(SDMQ25==SDMQ75,SDMQ25-0.02,SDMQ25)
        SDMQ25<-ifelse(SDMQ25==SDMQ50,SDMQ25-0.01,SDMQ25)
        SDMQ25<-ifelse(SDMQ25<0.0000001,0.0000001,SDMQ25)
        print(paste("SDMQ25 =",SDMQ25))
        SDMqs<-cutSDM

        if(is.na(SDMQ75)){
          SDMQ75<-0.75
          SDMQ50<-0.5
          SDMQ25<-0.25
          SDMqs<-KTO
          print(paste("default SDMQ = 0.25, 0.5, 0.75"))

        }

        values(SDMqs)<-ifelse(values(SDMqs)>Q75, 4, values(SDMqs))
        values(SDMqs)<-ifelse(values(SDMqs)>Q50 & values(SDMqs) <= Q75, 3, values(SDMqs))
        values(SDMqs)<-ifelse(values(SDMqs)>Q25 & values(SDMqs) <= Q50, 2, values(SDMqs))
        values(SDMqs)<-ifelse(values(SDMqs)>0.0000001 & values(SDMqs) <= Q25, 1, values(SDMqs))
        values(SDMqs)<-ifelse(values(SDMqs)<=0.0000001, NA, values(SDMqs))

        Sys.sleep(0.1)
        print(paste(Sys.time(), "4b. coverting to shapefile", sep=" "))

        SDMqsshp<-rasterToPolygons(SDMqs, dissolve=TRUE)
        SDMqsshp<-st_as_sf(SDMqsshp)
        names(SDMqsshp)[1]<-"VALUE"
        SDMqsshp$TYPE<-rep("Suitability Quantile Categories",nrow(SDMqsshp))
        SDMqsshp$SPECIES<-rep(gsub("_"," ",sp),nrow(SDMqsshp))
        SDMqsshp$MODUNIT<-rep(gsub("_"," ",mapsp),nrow(SDMqsshp))
        SDMqsshp$CODE_SSN<-rep(DatSum$SSN[n],nrow(SDMqsshp))
        SDMqsshp$DATETIME<-rep(date(),nrow(SDMqsshp))
        SDMqsshp$CAT<-rep(cat,nrow(SDMqsshp))
        st_write(SDMqsshp, paste(IndSppOutDir,"/",sp,"_",cat,"_SDMQs.shp",sep=""),delete_layer=TRUE)

        
        #########################################
        #categorized suitability:
        
        
        Sys.sleep(0.1)
        print(paste(Sys.time(), "5a. SDM equal categories", sep=" "))
        
        SDMcat<-cutSDM
        SDMQ75<-round(quantile(myvals,0.75, na.rm=TRUE),2)
        
        if(is.na(SDMQ75)){
          SDMcat<-KTO
        }
        
        # values(SDMcat)<-ifelse(values(SDMcat)<=0.0000001, NA, values(SDMcat))
        # myvals<-values(cutSDM)
        # SDMQ75<-round(quantile(myvals,0.75, na.rm=TRUE),2)
        # if(is.na(SDMQ75)){
        #   SDMcat<-KTOBuf
        # }
        # 
        values(SDMcat)<-ifelse(values(SDMcat)>0.75, 4, values(SDMcat))
        values(SDMcat)<-ifelse(values(SDMcat)>0.5 & values(SDMcat) <= 0.75, 3, values(SDMcat))
        values(SDMcat)<-ifelse(values(SDMcat)>0.25 & values(SDMcat) <= 0.5, 2, values(SDMcat))
        values(SDMcat)<-ifelse(values(SDMcat)>0.0000001 & values(SDMcat) <= 0.25, 1, values(SDMcat))
        values(SDMcat)<-ifelse(values(SDMcat)<=0.0000001, NA, values(SDMcat))
        
        Sys.sleep(0.1)
        print(paste(Sys.time(), "5b. coverting to shapefile", sep=" "))
        
        SDMcatshp<-rasterToPolygons(SDMcat, dissolve=TRUE)
        SDMcatshp<-st_as_sf(SDMcatshp)
        names(SDMcatshp)[1]<-"VALUE"
        SDMcatshp$TYPE<-rep("Suitability Value Categories",nrow(SDMcatshp))
        SDMcatshp$SPECIES<-rep(gsub("_"," ",sp),nrow(SDMcatshp))
        SDMcatshp$MODUNIT<-rep(gsub("_"," ",mapsp),nrow(SDMcatshp))
        SDMcatshp$CODE_SSN<-rep(DatSum$SSN[n],nrow(SDMcatshp))
        SDMcatshp$DATETIME<-rep(date(),nrow(SDMcatshp))
        SDMcatshp$CAT<-rep(cat,nrow(SDMcatshp))
        st_write(SDMcatshp, paste(IndSppOutDir,"/",sp,"_",cat,"_SDMcat.shp",sep=""),delete_layer=TRUE)
        
      }
    }
  }
}


######################################################
#delete the temp dir for this scriptL
unlink(paste(MyDir,"/RastoShape",sep=""), recursive = TRUE) 
######################################################

